double  interaccion(caja *bx)
{
	int ind;
	int ank, ank2, ankh;
	double virial=0.0;
	ank=6;
	ank2=2*ank;
	ankh=ank/2;
	double dist2,rqInv,dr[3], tsn, ccel3,potdm, ePoten;
	double conrep= 4.0; //Parte repulsiva del potencial Leonard-Jones
	double conref= 48.0; //Parte repulsiva de la fuerza del Leonard-Jones
	double conat= 4.0; //Parte atractiva del potencial Leonard-Jones
	double conatf= 24.0;// Parte atractiva de la fuerza del Leonard-Jones
	double rc = bx->rcut;
	double rc2=rc*rc;

	//Ponemos las fuerzas de las particulas de la caja a 0
	bx->initForces();

	//Correccion al cutoff del potencial
	double LJcut= pow(1.0/(rc*rc),ankh)*(pow(1.0/(rc*rc),ankh)*conrep-conat);

	ePoten=0.;
	int count =0;
	for(int i = 0;i<bx->npart-1;i++)
	{
		for(int j=i+1;j<bx->npart;j++)
		{

			//////////////
			///Calculamos la distancia entre i y j
			dist2= bx->distancia2(i,j);
			if(dist2<rc2)
			{
				count++;
				rqInv =1.0/dist2;//1/dist^2
				tsn=pow(rqInv,ankh);//1/dist^6
				//ccel3 = F/r=-(dPOTDM/dr)/r 
				ccel3 =rqInv*tsn*(conref*tsn-conatf);
				potdm = tsn*(conrep*tsn-conat);
				//Acualizamos la energia potencial
				ePoten+=potdm;//-LJcut;
				//Calculamos los componentes de las fuerzas
				for (int k=0;k<3;k++)
				{
					dr[k]= bx->dr[k];
					//Para la particula iii
					bx->part[i].force[k]+=ccel3*dr[k];
					//Para la particula jjj
					bx->part[j].force[k]-=ccel3*dr[k];
				}
				//Actualizamos el virial
				virial-=ccel3*(dr[0]*dr[0]-dr[1]*dr[1]-dr[2]*dr[2]);
			}
		}
	}
	//Introducimos la correccion del cutoff del potencial
	ePoten-=LJcut*count;//(bx->npart)*((bx->npart)-1)/2.0;
	bx->epoten=ePoten/double(bx->npart);
	bx->virial=  virial;
	//bx->virial=0.95*bx->vol;//-virial/3.0;
	return ePoten;

}
